{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%202%20-%20Creating%20and%20Simulating%20a%20simple%20arm%20Model.ipynb\">\n",
    "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "p60SkX5SZXDS"
   },
   "source": [
    "# 2. Creating and Simulating a Simple Arm Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YLK6vcNo1NHP"
   },
   "source": [
    "## 2.1. Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TKiEOjfi1LW1"
   },
   "source": [
    "**Purpose**\n",
    "\n",
    "The purpose of this tutorial is to introduce you to the creation of musculoskeletal models using [OpenSim](https://opensim.stanford.edu/) [[1]](https://doi.org/10.1109/TBME.2007.901024). In this tutorial you will:\n",
    "\n",
    "*   Create a simple arm model.\n",
    "*   Use the model to perform a forward simulation.\n",
    "*   Plot information obtained from the simulation.\n",
    "\n",
    "**Format**\n",
    "\n",
    "In this tutorial, you will create an arm model. The model is very simple and consists of two bones (humerus and radius), and a muscle to flex the arm (biceps).\n",
    "\n",
    "The first step of this tutorial will be to set up the environment, and then you will programmatically create the model by defining the necessary bodies, muscles, and joints. Finally, you will perform Forward Kinematics [[2]](https://mitpress.mit.edu/9780262044202/biomechanics-of-movement/) with the model to obtain kinematic data.\n",
    "\n",
    "**Local Installation**\n",
    "\n",
    "You may also run this notebook using a locally installed OpenSim `conda` package. Visit [the OpenSim python scripting Confluence page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python) for more information on how to install OpenSim in a local Python environment. Then, you may skip to Section 2.3 below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ADllcxbsjTr8"
   },
   "source": [
    "## 2.2. Setup `conda` and OpenSim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mmmTuEI4HQJG"
   },
   "source": [
    "First, set up the environment by executing the following cell (See [Tutorial 1: Introduction to OpenSim](https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%201%20-%20Introduction%20to%20OpenSim.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kdnS0S_xioqg"
   },
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vt5loGbdfKb8"
   },
   "source": [
    "Don't worry if after executing the previous cell you get an error saying that your session has failed. This happens because `condacolab` needs to restart the session for the changes to take effect. Therefore, you have to execute the previous cell before executing any other cell.\n",
    "\n",
    "Now, execute the following cell to install the OpenSim `conda` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6ZnEPfPgjup9"
   },
   "outputs": [],
   "source": [
    "!conda install opensim-org::opensim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Q85J8kOTj69D"
   },
   "source": [
    "## 2.3. Create the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "U_hlEyz5HT4E"
   },
   "source": [
    "First, import OpenSim and create an empty `Model` named \"simple_arm\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ChnLbJ0CkSH7"
   },
   "outputs": [],
   "source": [
    "import opensim as osim\n",
    "\n",
    "# Create empty model and assign to it the name \"simple_arm\".\n",
    "arm = osim.Model()\n",
    "arm.setName(\"simple_arm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ENIL146IlcUk"
   },
   "source": [
    "Next, create two bodies (the bones) using the `Body` class constructor function ([documentation](https://simtk.org/api_docs/opensim/api_docs/classOpenSim_1_1Body.html#ab1d061783d6e0b093e0794c19326e7c0)). This function takes 4 inputs:\n",
    "  * **name**: use \"humerus\" and \"radius\" for the names of the two bodies\n",
    "  * **mass**: set a mass of 1 $kg$ for both bodies\n",
    "  * **center of mass location**: set the location to be the body's origin (i.e., $(0,0,0)$) in both bodies (`osim.Vec3(0)`)\n",
    "  * **inertia**: use an inertia matrix with moments of inertia of 0.1, 0.5 and 0.1 $kg \\cdot m^2$, in the x, y and z directions, respectively, and with products of inertia of 0 (`osim.Inertia(0.1, 0.5, 0.1)`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JP4vs1OtlfOa"
   },
   "outputs": [],
   "source": [
    "# Define the humerus.\n",
    "humerus = osim.Body(\"humerus\",\n",
    "                    1.0,\n",
    "                    osim.Vec3(0),\n",
    "                    osim.Inertia(0.1, 0.5, 0.1))\n",
    "\n",
    "# Define the radius.\n",
    "radius = osim.Body(\"radius\",\n",
    "                   1.0,\n",
    "                   osim.Vec3(0),\n",
    "                   osim.Inertia(0.1, 0.5, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JZM8_weblhfc"
   },
   "source": [
    "The model needs pin joints, and we'll use the `PinJoint` class constructor function (which uses the same input as `Joint`, [documentation](https://simtk.org/api_docs/opensim/api_docs/classOpenSim_1_1Joint.html#a31854c22f626c982ad67fcff12a94961)) to create two pin joints. The first pin joint corresponds with the shoulder, and connects the top of the humerus with the ground (the humerus is the root of the model). The second pin joint corresponds with the elbow and connects the bottom of the humerus with the top of the radius. Assume each body is 1 $m$ long."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QxAX4pUrli0N"
   },
   "outputs": [],
   "source": [
    "# Shoulder pin joint from the ground to the humerus, so it is fixed and will not fall.\n",
    "shoulder = osim.PinJoint(\"shoulder\",\n",
    "                         arm.getGround(), # PhysicalFrame\n",
    "                         osim.Vec3(0),\n",
    "                         osim.Vec3(0),\n",
    "                         humerus, # PhysicalFrame\n",
    "                         osim.Vec3(0, 0.5, 0), # attach the joint 0.5 m up from the humerus origin\n",
    "                         osim.Vec3(0))\n",
    "\n",
    "# Elbow pin joint from the humerus to the radius.\n",
    "elbow = osim.PinJoint(\"elbow\",\n",
    "                      humerus, # PhysicalFrame\n",
    "                      osim.Vec3(0, -0.5, 0), # attach the joint 0.5 m down from the humerus origin\n",
    "                      osim.Vec3(0),\n",
    "                      radius, # PhysicalFrame\n",
    "                      osim.Vec3(0, 0.5, 0), # attach the joint 0.5 m up from the radius origin\n",
    "                      osim.Vec3(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MmIGq-lYlk2W"
   },
   "source": [
    "In order to flex the elbow, now you will add a muscle representing the biceps to the model. This muscle is an actuator that can be actively controlled to generate force. \n",
    "\n",
    "The muscle is created as a configurable muscle using the `Millard2012EquilibriumMuscle` class [[3]](https://doi.org/10.1115/1.4023390), and has a maximum isometric force of 100 $N$, optimal fiber length of 0.6 $m$, tendon slack length of 0.55 $m$, and a pennation angle of 0.0 $rad$. Then, two path points are defined with muscle's origin point in the humerus and insertion point in the radius.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BkNnZOJtllMu"
   },
   "outputs": [],
   "source": [
    "# Create and set the parameters for the biceps muscle.\n",
    "biceps = osim.Millard2012EquilibriumMuscle(\"biceps\",  # Muscle name\n",
    "                                           100.0,     # Max isometric force\n",
    "                                           0.6,       # Optimal fiber length\n",
    "                                           0.55,      # Tendon slack length\n",
    "                                           0.0)       # Pennation angle\n",
    "\n",
    "# Add path points to the humerus and radius. The allows the muscle to generate\n",
    "# forces on these two bodies.\n",
    "biceps.addNewPathPoint(\"origin\",\n",
    "                       humerus,\n",
    "                       osim.Vec3(0, 0.3, 0))\n",
    "biceps.addNewPathPoint(\"insertion\",\n",
    "                       radius,\n",
    "                       osim.Vec3(0, 0.2, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "w4dbDoA6loHI"
   },
   "source": [
    "In order to excite the muscle (activate the actuator), the model needs a controller. In the following cell you will create a `PrescribedController`, add the biceps as an actuator to the controller, and specify the excitation of the muscle. In this example, the excitation is time-varying and defined by a step function ([documentation](https://simtk.org/api_docs/opensim/api_docs/classOpenSim_1_1StepFunction.html#afdce40906fc61bc7dbcfa026d1f0ad8b)).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "IKcMjiBIlpY2"
   },
   "outputs": [],
   "source": [
    "# Create a controller, add the actuator (muscle) to it, and specify the excitation.\n",
    "brain = osim.PrescribedController()\n",
    "brain.addActuator(biceps)\n",
    "brain.prescribeControlForActuator(\"biceps\", osim.StepFunction(0.5, 3.0, 0.3, 1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9YB5snO5lruj"
   },
   "source": [
    "Now that every component of the model has been defined, let's add all of the elements to the arm model. Remember that the elements that you have defined are the bodies (bones) of the model, the joints of the model, an actuator (muscle) that generates force in the model, and a controller (brain) that specifies how the actuator is activated throughout the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "glSBsxAaltLg"
   },
   "outputs": [],
   "source": [
    "# Add all of the elements defined in previous cells to the model.\n",
    "arm.addBody(humerus)\n",
    "arm.addBody(radius)\n",
    "arm.addJoint(shoulder)\n",
    "arm.addJoint(elbow)\n",
    "arm.addForce(biceps)\n",
    "arm.addController(brain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wGKY1uv2l3vt"
   },
   "source": [
    "After building the model, you have to initialize the system in order to access and change the state of the model. \n",
    "\n",
    "We lock the shoulder joint so that it will not change its angle during the simulation. We set the elbow angle to an initial value of 90º, or `0.5 * π` $rad$. Finally, we update the state of the muscle, so that the force in the muscle and in the tendon are in equilibrium before we start the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "npMwjT5ol3_V"
   },
   "outputs": [],
   "source": [
    "# Initialize model.\n",
    "state = arm.initSystem()\n",
    "\n",
    "# Fix the shoulder at its default angle and begin with the elbow flexed 90º degrees.\n",
    "shoulder.getCoordinate().setLocked(state, True)\n",
    "elbow.getCoordinate().setValue(state, 0.5 * osim.SimTK_PI)\n",
    "arm.equilibrateMuscles(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MqQqxE20l-Yd"
   },
   "source": [
    "Now, you can simulate the model. First, let's create a `Manager` and set the starting time to the state. The `Manager` is an object that manages the execution of the simulation. Using the manager, you can then perform a \"forward simulation\" [[2]](https://mitpress.mit.edu/9780262044202/biomechanics-of-movement/), where the manager integrates the equations of motion for the specified model forward in time, given the current state and a final time. This will return a new state with the results of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vBpKDPgsl_Zh"
   },
   "outputs": [],
   "source": [
    "# Initialize the system and set initial time.\n",
    "state.setTime(0)\n",
    "\n",
    "# Get and initialize manager with the state.\n",
    "manager = osim.Manager(arm)\n",
    "manager.initialize(state)\n",
    "\n",
    "# Integrate the model until the ending time.\n",
    "state = manager.integrate(10.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XEumXmZyWQ2b"
   },
   "source": [
    "After the simulation is finished, you can extract the table of states. Let's also print the column labels, so we have the paths to the states of the values that we are interested in plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mbLScbezWa2W"
   },
   "outputs": [],
   "source": [
    "# Extract table of states.\n",
    "statesTable = manager.getStatesTable()\n",
    "\n",
    "# Print labels of columns in the table.\n",
    "print(statesTable.getColumnLabels())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KqD7waFqRo_W"
   },
   "source": [
    "We want to plot the shoulder and elbow positions and speed, as well as the biceps activation and fiber length over time. Let's extract the desired data from the table of states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HNWVlK-VRp2C"
   },
   "outputs": [],
   "source": [
    "# Extract time values for the x-axis.\n",
    "times = statesTable.getIndependentColumn()\n",
    "\n",
    "# Extract shoulder data.\n",
    "shoulder_pos = statesTable.getDependentColumn('/jointset/shoulder/shoulder_coord_0/value')\n",
    "shoulder_speed = statesTable.getDependentColumn('/jointset/shoulder/shoulder_coord_0/speed')\n",
    "\n",
    "# Extract elbow data.\n",
    "elbow_pos = statesTable.getDependentColumn('/jointset/elbow/elbow_coord_0/value')\n",
    "elbow_speed = statesTable.getDependentColumn('/jointset/elbow/elbow_coord_0/speed')\n",
    "\n",
    "# Extract biceps data.\n",
    "biceps_activation = statesTable.getDependentColumn('/forceset/biceps/activation')\n",
    "biceps_fiber_length = statesTable.getDependentColumn('/forceset/biceps/fiber_length')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PAJ0FRKHardS"
   },
   "source": [
    "## 2.4. Plotting data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "D-iFifvcVR5U"
   },
   "source": [
    "In this section, you will plot the results of the simulation using a well known python library: [Matplotlib](https://matplotlib.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kX6MdpfEa3aS"
   },
   "source": [
    "### 2.4.1. Matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-Vu58PadWMQh"
   },
   "source": [
    "The following cell plots the desired columns (`shoulder_pos`, `shoulder_speed`, `elbow_pos`, `elbow_speed`, `biceps_act` and `biceps_fib_len`) using `matplotlib`. Each column is plotted in a different subplot.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HuquEKo-a4_E"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Create six subplots, with 2 rows and 3 columns.\n",
    "fig, axs = plt.subplots(2, 3, figsize=(14, 8))\n",
    "fig.suptitle('Arm model data representation')\n",
    "\n",
    "# Plot shoulder position and speed in column 1.\n",
    "axs[0, 0].plot(times, shoulder_pos.to_numpy())\n",
    "axs[0, 0].set_title('Shoulder pos.')\n",
    "axs[0, 0].set_xlabel('Time')\n",
    "axs[0, 0].set_ylabel('Position')\n",
    "axs[0, 0].grid()\n",
    "axs[1, 0].plot(times, shoulder_speed.to_numpy())\n",
    "axs[1, 0].set_title('Shoulder speed')\n",
    "axs[1, 0].set_xlabel('Time')\n",
    "axs[1, 0].set_ylabel('Speed')\n",
    "axs[1, 0].grid()\n",
    "\n",
    "# Plot elbow position and speed in column 2.\n",
    "axs[0, 1].plot(times, elbow_pos.to_numpy())\n",
    "axs[0, 1].set_title('Elbow pos.')\n",
    "axs[0, 1].set_xlabel('Time')\n",
    "axs[0, 1].set_ylabel('Position')\n",
    "axs[0, 1].grid()\n",
    "axs[1, 1].plot(times, elbow_speed.to_numpy())\n",
    "axs[1, 1].set_title('Elbow speed')\n",
    "axs[1, 1].set_xlabel('Time')\n",
    "axs[1, 1].set_ylabel('Speed')\n",
    "axs[1, 1].grid()\n",
    "\n",
    "# Plot biceps activation and length in column 3.\n",
    "axs[0, 2].plot(times, biceps_activation.to_numpy())\n",
    "axs[0, 2].set_title('Biceps act.')\n",
    "axs[0, 2].set_xlabel('Time')\n",
    "axs[0, 2].set_ylabel('Activation')\n",
    "axs[0, 2].grid()\n",
    "axs[1, 2].plot(times, biceps_fiber_length.to_numpy())\n",
    "axs[1, 2].set_title('Biceps fiber length')\n",
    "axs[1, 2].set_xlabel('Time')\n",
    "axs[1, 2].set_ylabel('Fiber length')\n",
    "axs[1, 2].grid()\n",
    "\n",
    "# Set the spacing between subplots\n",
    "plt.subplots_adjust(wspace=0.5, hspace=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "i7rhPT2id6gg"
   },
   "source": [
    "## 2.5. Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1XK9UgItWQMF"
   },
   "source": [
    "In this tutorial you have created a simple arm model defined by two bones (bodies), two joints, and a muscle (actuator). Then, you have simulated the model and plotted the calculated values.\n",
    "\n",
    "To create more complex models graphically and visualize them, check out [OpenSim Creator](https://opensimcreator.com/). OpenSim Creator is freely available, open-source software that enables users to visually build, edit, and simulate OpenSim models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lLWgd8Noun6Y"
   },
   "source": [
    "## 2.6. Useful Links\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m8M3dlT5xQcD"
   },
   "source": [
    "> **OpenSim Website:** https://opensim.stanford.edu/\n",
    ">\n",
    "> **OpenSim Python Scripting:** https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python\n",
    ">\n",
    "> **OpenSim API Documentation:** https://simtk.org/api_docs/opensim/api_docs/\n",
    "> \n",
    "> **OpenSim Creator Website:** https://opensimcreator.com/\n",
    "> \n",
    "> **SimTK Website:** https://simtk.org/projects/opensim\n",
    "> \n",
    "> **Biomechanics of Movement Course Videos:** https://www.youtube.com/channel/UCDNGy0KKNLQ-ztcL5h2Z6zA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rFR2Rf9tk_9A"
   },
   "source": [
    "## 2.7 Acknowledgments\n",
    "\n",
    "Thanks to [OpenSimColab](https://simtk.org/projects/opencolab) project [[4]](https://doi.org/10.1080/10255842.2022.2104607) for creating the first OpenSim Conda package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dxMjdukAnH2N"
   },
   "source": [
    "## 2.8. References\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WX2A_iaYulV1"
   },
   "source": [
    "> [1] Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., Guendelman, E., & Thelen, D. G. (2007). **OpenSim: open-source software to create and analyze dynamic simulations of movement.** *IEEE Transactions on Bio-Medical Engineering*, 54(11), 1940–1950. https://doi.org/10.1109/TBME.2007.901024\n",
    ">\n",
    "> [2] Uchida, T. K., Delp, S. L., & Delp, D. (2021). **Biomechanics of movement: The science of sports, robotics, and rehabilitation.** *MIT Press*. https://mitpress.mit.edu/9780262044202/biomechanics-of-movement/\n",
    ">\n",
    "> [3] Millard, M., Uchida, T., Seth, A., & Delp, S. L. (2013). **Flexing computational muscle: Modeling and simulation of musculotendon dynamics.** *Journal of Biomechanical Engineering*, 135(2). https://doi.org/10.1115/1.4023390\n",
    ">\n",
    "> [4] Mokhtarzadeh, H., Jiang, F., Zhao, S., & Malekipour, F. (2022). **OpenColab project: OpenSim in Google colaboratory to explore biomechanics on the web.** *Computer Methods in Biomechanics and Biomedical Engineering*, 1–9. https://doi.org/10.1080/10255842.2022.2104607"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "opensim",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
